function v_mu = mos_variance(Z)
% Determines the variance of the mean opinion score, given a matrix of scores.
% Unknown scores are represented by NaN.
%
% Z: a N-by-M matrix of scores, where the rows are subjects and columns are sentences
% We assume that
%
% 	Z_ij = mu + x_i + y_j + eps_ij, where
%
% mu is the mean opinion score (given by nanmean(Z(:)))
% x_i ~ N(0, sigma_w^2), with sigma_w^2 modeling worker variation
% y_j ~ N(0, sigma_s^2), with sigma_y^2 modeling sentence variation
% eps_ij ~ N(0, sigma_u^2), with sigma_u^2 modeling worker uncertainty
%
% The returned value v_mu is Var[mu].

[N, M] = size(Z);

W  = ~isnan(Z);
Mi = sum(W,1); Mi = Mi(:);
Nj = sum(W,2); Nj = Nj(:);
T  = sum(W(:));

v_su  = vertical_var(Z');	% v_su  = v_s + v_u
v_wu  = vertical_var(Z);	% v_wu  = v_w + v_u
v_swu = nanvar(Z(:));		% v_swu = v_s + v_w + v_u
if ~isnan(v_su) && ~isnan(v_wu)
	v = [1 0 1 ; 0 1 1 ; 1 1 1] \ [ v_su ; v_wu ; v_swu ];		% v = [v_s ; v_w ; v_u]
	v = max(v,0);

	v_s = v(1);
	v_w = v(2);
	v_u = v(3);

	v_mu = v_s * sum(Mi.^2)/T^2 + v_w * sum(Nj.^2)/T^2 + v_u/T;
elseif isnan(v_su) && ~isnan(v_wu)
	v = [0 1 ; 1 1] \ [ v_wu ; v_swu ];		% v = [v_s ; v_wu]
	v = max(v,0);

	v_s  = v(1);
	v_wu = v(2);

	v_mu = v_s * sum(Mi.^2)/T^2 + v_wu/T;
elseif ~isnan(v_su) && isnan(v_wu)
	v = [0 1 ; 1 1] \ [ v_su ; v_swu ];		% v = [v_w ; v_su]
	v = max(v,0);

	v_w  = v(1);
	v_su = v(2);

	v_mu = v_w * sum(Nj.^2)/T^2 + v_su/T;
else
	assert(isnan(v_su) && isnan(v_wu));
	v_mu = v_swu/T;
end

%%%%

function v = vertical_var(Z)
[rows, cols] = size(Z);
v = [];
for i = 1:cols
	if nancount(Z(:,i)) >= 2
		v = [v ; nanvar(Z(:,i))];
	end
end
v = mean(v);

%%%%

function n = nancount(vec)
n = sum(~isnan(vec));

